%% Cholesky and Proxy SVARs in Levels with Bias Correction
%  Mian Sufi Verner, 2017/03/06
%% Figure 1, Figure 4, Figure A1, Figure A2
% (1) Define and estimate reduced form panel VAR in levels with country FE
% (2) Bootstrap to estimate bias and bias correction term in VAR betas
% (3) Bootstrap to construct Cholesky IRF confidence intervals (bias corrected)
% (4) Bootstrap to construct Proxy IRF confidence intervals (bias corrected)
%==========================================================================

clc; clear all; close all;
rng default
set(0,'defaulttextfontsize',17);
set(0,'defaultaxesfontsize',17);
set(0,'defaultlinelinewidth',2);
ifig=1;



%% SVAR.(HHD(t)/Y(t-1),FD(t)/Y(t-1),lnY(t)). 
%  Proxy SVAR Instrument = dummy for low mortgage spread 

%% (0) Prepare Data 
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

[STATA]=csvread('VAR3_levels_fixedeffects_IVprctiles.csv');
i_m = 3; I_Y = 7:1:9; I_X = 10:1:size(STATA,2);                            
I_mNaN     = find(STATA(:,i_m)==999999);                                   % (missing coded 999999)
VAR.I_m    = find(STATA(:,i_m)~=999999);
STATA(I_mNaN,i_m)=nan;
res_instrument = 1;                                                        % use residual instrument (residual on instrument on X)
VAR.m         = STATA(:,i_m);
cc            = STATA(:,1);                                                
[VAR.Ti,VAR.ii]= histc(cc,unique(cc));                                     
VAR.Year      = STATA(:,2);
VAR.Y         = STATA(:,I_Y); 
VAR.X         = STATA(:,I_X);                                              
VAR.p         = 5;                                                         % number of lags
VAR.irhor     = 11;                                                       
VAR.h         = 0:1:(VAR.irhor-1);
VAR.N         = max(VAR.ii);                                               
VAR.Perm      = [0 0 1; 0 1 0 ; 1 0 0];                                    % Cholesky ordering = lnY, FD, HHD
VAR.normalize = 1;
clevel        = [95 ; 90];

% orthogonalize instrument with respect to X
if res_instrument == 1;
    Xrtmp     = [VAR.X(VAR.I_m,:) ones(length(VAR.X(VAR.I_m)),1)]; 
    gamma     = Xrtmp\VAR.m(VAR.I_m,1);  
    VAR.m(VAR.I_m,1) = VAR.m(VAR.I_m,1) - Xrtmp*gamma;
end



%% (1) Estimate RF VAR
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
VAR      = rfvar(VAR);

%% (2) Bias correct RF VAR
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

% step 0: create D (panel dummy variable matrix), demeaned Y, demeaned X
D = zeros(VAR.T,VAR.N);
t  = 1;
for i=1:VAR.N;
    D(t:(t+VAR.Ti(i)-1),i) = ones(VAR.Ti(i),1);
                         t = t + VAR.Ti(i);
end
% demean Y and X
b_means    = (D'*D)\(D'*[VAR.Y VAR.X(:,1:(VAR.Ny*VAR.p))]);
VAR.Ydm    = VAR.Y - D*b_means(:,1:VAR.Ny);
VAR.Xdm    = VAR.X(:,1:(VAR.Ny*VAR.p)) - D*b_means(:,(VAR.Ny+1):end);
VAR.bettil = VAR.Xdm \ VAR.Ydm;
VAR.D      = D;

% step 1: define inputs to bias correction function
tol          = 0.0015; 
nboot        = 2000;
delta_update = .5;
VAR.bet0 = VAR.bet;                                                        % initial guess for bet_bc is OLS beta

% step 2: run bias correction function. 
VARbc = doVAR_BiasCorrected(VAR,nboot,clevel,tol,delta_update);

% step 3: extract bias and bias corrected beta from VARbc 
VAR.bc_diff = VAR.bet(1:(VAR.Ny*VAR.p),:) - VARbc.bet(1:(VAR.Ny*VAR.p),:); % bias 
VAR.bet_bc  = [VARbc.bet ; VAR.bet( (VAR.Ny*VAR.p+1:end),:)];              % bias corrected beta

% adjust fixed effects 
    uhat     = VAR.Y - VAR.X(:,1:VAR.Ny*VAR.p)*VAR.bet_bc(1:VAR.Ny*VAR.p,:);
    D1       = [D(:,1:(VAR.N-1)) ones(VAR.T,1)];
    alpha_bc = (D1'*D1)\(D1'*uhat);
VAR.bet_bc(VAR.Ny*VAR.p+1:end,:) = alpha_bc;
VAR.res_bc   = VAR.Y - [VAR.X ones(VAR.T,1)]*VAR.bet_bc;
VAR.Sigma_bc = VAR.res_bc'*VAR.res_bc / (VAR.T - VAR.Nx);                  % bias corrected Sigma
VAR.res      = VAR.res_bc;                                                 % update residuals

% residuals need to be rescaled to
% account for fact that estimated error terms have lower variance than pop
% error terms (de Vos, Everaert, Ruyseen eqn (8))
VAR.res     = VAR.res * sqrt(VAR.T/(VAR.T-VAR.Nx)); 
VAR.Sigma   = VAR.res'*VAR.res/(VAR.T - VAR.Nx);
 

%% (3) Choleksy VAR with bias corrected beta and Sigma
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

% uncorrected Cholesky IRF
h   = VAR.h;
VAR = irfvar(VAR);

% bias corrected Cholesky IRF
VARbc.Ny    = VAR.Ny; VARbc.p=VAR.p; VARbc.Perm = VAR.Perm;
VARbc.irhor = VAR.irhor; VARbc.normalize = VAR.normalize;
VARbc.bet   = VAR.bet_bc;
VARbc.Sigma = VAR.Sigma_bc;
VAR_Chol_bc  = irfvar(VARbc);

% Figure A1
figure(ifig);
subplot(1,3,1); plot(h,VAR_Chol_bc.cholIRF(:,1,1),'o-k',h,VAR.cholIRF(:,1,1),'s--k'); title('d^{HH} \rightarrow d^{HH}'); legend({'Bias corrected','Not bias corrected'});
subplot(1,3,2); plot(h,VAR_Chol_bc.cholIRF(:,3,1),'o-k',h,VAR.cholIRF(:,3,1),'s--k'); title('d^{HH} \rightarrow y'); axis([0,10,-.5,.4]);
subplot(1,3,3); plot(h,VAR_Chol_bc.cholIRF(:,3,2),'o-k',h,VAR.cholIRF(:,3,2),'s--k'); title('d^{F} \rightarrow y'); axis([0,10,-.5,.4]);
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 .55]);
ifig=ifig+1;




%% (4) Choleksy VAR with bootstrapped CIs, bias corrected
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

VAR_Chol_bs_bc = doCholVARbootstrap_BiasCorrected_CS(VAR,nboot,clevel);



% Figure 1
figure(ifig); 
subplot(1,3,1); plot(h,VAR_Chol_bc.cholIRF(:,1,1),'o-k',h,VAR_Chol_bs_bc.CholirsH1(:,1,1),'--k',h,VAR_Chol_bs_bc.CholirsL1(:,1,1),'--k'); title('d^{HH} \rightarrow d^{HH}');
legend({'Impulse response','95% CI'},'Location','SouthWest');  
subplot(1,3,2); plot(h,VAR_Chol_bc.cholIRF(:,3,1),'o-k',h,VAR_Chol_bs_bc.CholirsH1(:,3,1),'--k',h,VAR_Chol_bs_bc.CholirsL1(:,3,1),'--k'); title('d^{HH} \rightarrow y'); xlabel('Years after impulse');
subplot(1,3,3); plot(h,VAR_Chol_bc.cholIRF(:,3,2),'o-k',h,VAR_Chol_bs_bc.CholirsH1(:,3,2),'--k',h,VAR_Chol_bs_bc.CholirsL1(:,3,2),'--k'); title('d^F \rightarrow y');
axis([0,10,-.5,.4]);
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 .55]);
ifig=ifig+1;


% Figuge A2: All 9 IRFS
figure(ifig); 
subplot(3,3,9); plot(h,VAR_Chol_bc.cholIRF(:,1,1),'o-k',h,VAR_Chol_bs_bc.CholirsH1(:,1,1),'--k',h,VAR_Chol_bs_bc.CholirsL1(:,1,1),'--k'); title('d^{HH} \rightarrow d^{HH}');
axis([0,10,0,2]);
subplot(3,3,8); plot(h,VAR_Chol_bc.cholIRF(:,2,1),'o-k',h,VAR_Chol_bs_bc.CholirsH1(:,2,1),'--k',h,VAR_Chol_bs_bc.CholirsL1(:,2,1),'--k'); title('d^{HH} \rightarrow d^F');
subplot(3,3,7); plot(h,VAR_Chol_bc.cholIRF(:,3,1),'o-k',h,VAR_Chol_bs_bc.CholirsH1(:,3,1),'--k',h,VAR_Chol_bs_bc.CholirsL1(:,3,1),'--k'); title('d^{HH} \rightarrow y');
axis([0,10,-.6,.6]);
subplot(3,3,6); plot(h,VAR_Chol_bc.cholIRF(:,1,2),'o-k',h,VAR_Chol_bs_bc.CholirsH1(:,1,2),'--k',h,VAR_Chol_bs_bc.CholirsL1(:,1,2),'--k'); title('d^F \rightarrow d^{HH}');
axis([0,10,-.3,.3]);
subplot(3,3,5); plot(h,VAR_Chol_bc.cholIRF(:,2,2),'o-k',h,VAR_Chol_bs_bc.CholirsH1(:,2,2),'--k',h,VAR_Chol_bs_bc.CholirsL1(:,2,2),'--k'); title('d^F \rightarrow d^F');
subplot(3,3,4); plot(h,VAR_Chol_bc.cholIRF(:,3,2),'o-k',h,VAR_Chol_bs_bc.CholirsH1(:,3,2),'--k',h,VAR_Chol_bs_bc.CholirsL1(:,3,2),'--k'); title('d^F \rightarrow y');
axis([0,10,-.6,.6]);
subplot(3,3,3); plot(h,VAR_Chol_bc.cholIRF(:,1,3),'o-k',h,VAR_Chol_bs_bc.CholirsH1(:,1,3),'--k',h,VAR_Chol_bs_bc.CholirsL1(:,1,3),'--k'); title('y \rightarrow d^{HH}');
subplot(3,3,2); plot(h,VAR_Chol_bc.cholIRF(:,2,3),'o-k',h,VAR_Chol_bs_bc.CholirsH1(:,2,3),'--k',h,VAR_Chol_bs_bc.CholirsL1(:,2,3),'--k'); title('y \rightarrow d^F');
subplot(3,3,1); plot(h,VAR_Chol_bc.cholIRF(:,3,3),'o-k',h,VAR_Chol_bs_bc.CholirsH1(:,3,3),'--k',h,VAR_Chol_bs_bc.CholirsL1(:,3,3),'--k'); title('y \rightarrow y');
legend({'Impulse response','95% CI'},'Location','SouthWest'); 
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 .75 1]);
ifig=ifig+1;



%% (5) Proxy VAR with bias corrected beta and Sigma
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

% not bias corrected
VAR_notbiascorr = doProxySVAR_single_EV(VAR);
% bias corrected
VAR             = doProxySVAR_biascorrected_single_EV(VAR);

figure(ifig);
subplot(2,2,1); plot(VAR.h,VAR.irs(:,1),'o-k',VAR.h,VAR_notbiascorr.irs(:,1),'s--k'); title('HHD_t/Y_{t-1} \rightarrow HHD_t/Y_{t-1}'); legend({'Bias corrected','Not bias corrected'});
subplot(2,2,2); plot(VAR.h,VAR.irs(:,2),'o-k',VAR.h,VAR_notbiascorr.irs(:,2),'s--k'); title('HHD_t/Y_{t-1} \rightarrow FD_t/Y_{t-1}');
subplot(2,2,3); plot(VAR.h,VAR.irs(:,3),'o-k',VAR.h,VAR_notbiascorr.irs(:,3),'s--k'); title('HHD_t/Y_{t-1} \rightarrow ln(Y_t)');
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 .7 1]);
ifig=ifig+1;


%% (6) Proxy VAR with bootstrapped CIs, bias corrected
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

nboot= 2000;
VAR.n = VAR.Ny;
% Get bootstrapped (bs), bias corrected (bc) impulse response function CIs
VAR_bs_bc = doProxySVARbootstrap_single_biascorrected_MSV_CS(VAR,nboot,clevel); 

% Figure 4: IRF with bootstrapped CIs: HHD->y, HHD->FD, HHD->HHD
figure(ifig);
subplot(1,3,2); plot(h,VAR.irs(:,3),'o-k',h,VAR_bs_bc.irsH1(:,3),'--k',h,VAR_bs_bc.irsL1(:,3),'--k'); title('d^{HH} \rightarrow y'); xlabel('Years after impulse');
axis([0, 10, -0.6, 0.6]);
subplot(1,3,3); plot(h,VAR.irs(:,2),'o-k',h,VAR_bs_bc.irsH1(:,2),'--k',h,VAR_bs_bc.irsL1(:,2),'--k'); title('d^{HH} \rightarrow d^{F}'); 
subplot(1,3,1); plot(h,VAR.irs(:,1),'o-k',h,VAR_bs_bc.irsH1(:,1),'--k',h,VAR_bs_bc.irsL1(:,1),'--k'); title('d^{HH} \rightarrow d^{HH}'); legend({'Impulse response','95% CI'},'Location','SouthWest');
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 .5]);
ifig=ifig+1;



